! to generate few group cross sections for initial MC guess
MODULE asmcal
    implicit none

    public init_asmcal, sync_asmarray, restore_core, setupasm_nxt, resetbuf_asmcal, homg_asm, pprec
    private findFASM

    type geoset
        sequence 
        real(8), pointer :: hx, hy
        integer, pointer :: cplst(:,:)
        real(8), pointer :: hacx(:), hacy(:), f2cx(:), f2cy(:)
        real(8) :: albedo(4)
        integer :: ncmx, ncmy, nfmx, nfmy
        integer, dimension(4) :: ibndry
    endtype

    integer, public :: ngcd
    
    
    integer, private, parameter :: max_ctinasm = 10
    type keyset
        sequence
        integer :: key
        integer :: value
    endtype
  
    
    integer :: tallymask
    integer, parameter :: tally_pinpower=1
    integer, parameter :: tally_cmfd=2
    integer, parameter :: tally_gcmfd=4
    integer, parameter :: tally_grpconst=8
    
    
!    type asmtallyset
!        sequence
!        ! i, j, g, thrd
!        real(8), pointer, dimension(:,:,:,:) :: phiv, rtot, rnuf, rabsr
!        ! multigroup -> few group
!        real(8), dimension(:) :: g2g
!    endtype
    
    real(8), pointer, private, dimension(:,:,:,:) :: asmphiv, asmtot, asmnuf, asmabs
    real(8), pointer, private, dimension(:,:,:,:,:) :: asmsct
    integer, private :: iasmtype=0
    
    type(keyset), pointer, dimension(:,:) :: xsec_match
    integer, pointer, public, dimension(:,:,:) :: xsec_map
    type(geoset), private :: asmgeo, coregeo, reflgeo
    
    
    type stxsec
        real(8), pointer, dimension(:,:) :: tot, nuf, abs, chi, rmv
        real(8), pointer, dimension(:,:,:) :: sct
        real(8), pointer, dimension(:,:) :: phiv
    endtype
    integer, private, pointer, dimension(:) :: maxct, maxctacc
    integer, private :: nmaxct
    integer, private ::ntot_ct=0
    logical, pointer, dimension(:) :: isfuelpin, isfuelasm
    
    type(stxsec) :: xsec_set
    real(8), private, pointer, dimension(:,:,:) :: shapef
contains 

subroutine init_asmcal
    include "pspec.FI"
    include "xsec.FI"
    integer :: iasm, i, j, k, maxctloc, ct
    
    
    ! temporarily setsk
    ngcd = gtog(ng)

    allocate(isfuelpin(ntypepin))
    allocate(isfuelasm(0:ntypeasm))
    isfuelpin = .false.
    isfuelasm = .false.

    
    allocate(asmphiv(ngcd,npx*2,npy, nthrd))
    allocate(asmtot(ngcd,npx*2,npy, nthrd))
    allocate(asmnuf(ngcd,npx*2,npy, nthrd))
    allocate(asmabs(ngcd,npx*2,npy, nthrd))
    allocate(asmsct(ngcd,ngcd,npx*2,npy, nthrd))
    call resetbuf_asmcal()

    allocate(xsec_match(max_ctinasm,ntypeasm))
    allocate(xsec_map(npx, npy, 0:ntypeasm))
    xsec_map=0
    allocate(maxct(ntypeasm))
    allocate(maxctacc(0:ntypeasm))
    maxct=0
    maxctacc=0
    
    
    if (homtype .eq. HT_PTASM) then
        ! 1. same cross section but different asm -> differnt xsec type
        ! 2. different cross section in an asm -> different 
        ! 3. therefore, existance search for the same asm
        ! setting xsec map
        do iasm=1, ntypeasm
            maxctloc=0
            do j=1, npy
                do i=1, npx
                    ct=assembly(i,j,iasm)
                    do k=1, maxctloc
                        if (ct .eq. xsec_match(k,iasm).key) exit
                    enddo
                    if (k>maxctloc) then
                        ntot_ct=ntot_ct+1
                        maxctloc=maxctloc+1
                        xsec_match(k,iasm).key=ct
                        xsec_match(k,iasm).value=ntot_ct
                        xsec_map(i,j,iasm)=ntot_ct
                    else
                        xsec_map(i,j,iasm)=xsec_match(k,iasm).value
                    endif
                enddo  
            enddo
            maxct(iasm)=maxctloc
            maxctacc(iasm)=maxctacc(iasm-1)+maxctloc
        enddo
        
        nmaxct=sum(maxct)
    elseif (homtype .eq. HT_ASMHM) then
        do iasm=1, ntypeasm
            ct=iasm
            ntot_ct=ntot_ct+1
            do j=1, npy
                do i=1, npx
                    xsec_match(1,iasm).key=ct
                    xsec_match(1,iasm).value=ntot_ct
                    xsec_map(i,j,iasm)=ntot_ct
                enddo
            enddo
            maxct(iasm)=1
            maxctacc(iasm)=maxctacc(iasm-1)+1
        enddo
        nmaxct=ntot_ct
    endif
    
    
!    allocate(xsec_set(nmaxct))
    allocate(xsec_set.tot(ngcd, nmaxct))
    allocate(xsec_set.nuf(ngcd, nmaxct))
    allocate(xsec_set.rmv(ngcd, nmaxct))
    allocate(xsec_set.abs(ngcd, nmaxct))
    allocate(xsec_set.chi(ngcd, nmaxct))
    allocate(xsec_set.sct(ngcd,ngcd, nmaxct))
    allocate(xsec_set.phiv(ngcd, nmaxct))
    
    if (homtype .eq. HT_ASMHM) then
        allocate(shapef(npx, npy, nmaxct))
        shapef=0
    endif

    xsec_set.tot=0
    xsec_set.nuf=0
    xsec_set.abs=0
    xsec_set.chi=0
    xsec_set.sct=0
    xsec_set.phiv=0
    
    ! save core input parameters
    coregeo.ncmx=ncmx
    coregeo.ncmy=ncmy
    coregeo.nfmx=nfmx
    coregeo.nfmy=nfmy
    coregeo.hacx=>hacx
    coregeo.hacy=>hacy
    coregeo.albedo=albedo
    coregeo.cplst=>cplst
    coregeo.ibndry=ibndry

    call findFASM()
endsubroutine

subroutine sync_asmarray(ithrd, phiv, tot, nuf, abs, sct)
    integer, intent(in) :: ithrd
    real(8), intent(inout), pointer, dimension(:,:,:) :: phiv, tot, nuf, abs
    real(8), intent(inout), pointer, dimension(:,:,:,:) :: sct
    
    phiv=>asmphiv(:,:,:,ithrd)
    tot=>asmtot(:,:,:,ithrd)
    nuf=>asmnuf(:,:,:,ithrd)
    abs=>asmabs(:,:,:,ithrd)
    sct=>asmsct(:,:,:,:,ithrd)
endsubroutine


subroutine resetbuf_asmcal
    include "pspec.FI"
    integer :: ithrd

    do ithrd=1, nthrd
        asmphiv(:,:,:,ithrd)=0
        asmtot(:,:,:,ithrd)=0
        asmnuf(:,:,:,ithrd)=0
        asmabs(:,:,:,ithrd)=0
        asmsct(:,:,:,:,ithrd)=0
    enddo
endsubroutine

!subroutine core2asm
!    implicit none
!    include "arrays.FI"
!    include "pspec.FI"
!    
!    coregeo.ncmx=ncmx
!    coregeo.ncmy=ncmy
!    coregeo.nfmx=nfmx
!    coregeo.nfmy=nfmy
!    coregeo.hacx=>hacx
!    coregeo.hacy=>hacy
!    coregeo.albedo=albedo
!    coregeo.cplst=>cplst
!    
!    asmgeo.ncmx=npx
!    asmgeo.ncmy=npx
!    allocate(asmgeo.hacx(0:asmgeo.ncmx))
!    allocate(asmgeo.hacy(0:asmgeo.ncmy))
!    allocate(asmgeo.cplst(asmgeo.ncmx,asmgeo.ncmy))
!    asmgeo.hacx=hacx(0:asmgeo.ncmx)
!    asmgeo.hacy=hacy(0:asmgeo.ncmy)
!    asmgeo.nfmx=asmgeo.ncmx
!    asmgeo.nfmy=asmgeo.ncmy
!    asmgeo.albedo=0.0_8
!    asmgeo.cplst(:,:)=assembly(:,:,iasmtype)
!    
!    ncmx=asmgeo.ncmx
!    ncmy=asmgeo.ncmy
!    nfmx=asmgeo.nfmx
!    nfmy=asmgeo.nfmy
!    albedo=asmgeo.albedo
!    cplst=>asmgeo.cplst
!    hacx=>asmgeo.hacx
!    hacy=>asmgeo.hacy
!
!    !allocate((0:nfmx))
!    
!endsubroutine    


subroutine homg_asm()
    include "pspec.FI"
    include "xsec.FI"
    integer :: ct
    integer :: i, j, k, g, gp, maxctloc, base, basei
    integer :: is=1
    type(stxsec), pointer :: xs
    character(256) :: sfile
    integer :: xt
    real(8) :: rsumpsi
    real(8), pointer, dimension(:,:) :: asmpsi
    

    
    if (not(isfuelasm(iasmtype))) then
        basei=npx
        ! is=npx/2
        is=1
    else
        basei=0
        is=1
        if (homtype .eq. HT_ASMHM) asmpsi=>shapef(:,:,iasmtype)
    endif
        
    do j=1, npy
        do i=is, npx
            ct=xsec_map(i,j, iasmtype)
            do g=1, ngcd
                xsec_set.phiv(g,ct)=xsec_set.phiv(g,ct)+sum(asmphiv(g,basei+i,j,:)) 
                xsec_set.tot(g,ct)=xsec_set.tot(g,ct)+sum(asmtot(g,basei+i,j,:)) 
                xsec_set.nuf(g,ct)=xsec_set.nuf(g,ct)+sum(asmnuf(g,basei+i,j,:)) 
                xsec_set.abs(g,ct)=xsec_set.abs(g,ct)+sum(asmabs(g,basei+i,j,:)) 
                do gp=1, ngcd
                    xsec_set.sct(gp,g,ct)=xsec_set.sct(gp,g,ct)+sum(asmsct(gp,g,basei+i,j,:)) 
                enddo
            enddo
            if (homtype .eq. HT_ASMHM .and. isfuelasm(iasmtype)) then
                asmpsi(i,j)=sum(asmnuf(:,basei+i,j,:))
            endif
        enddo
    enddo 
    
    base=maxctacc(iasmtype-1)
    do k=1, maxct(iasmtype)
        ct=base+k
        do g=1, ngcd
            xsec_set.tot(g,ct)=xsec_set.tot(g,ct)/xsec_set.phiv(g,ct)
            xsec_set.nuf(g,ct)=xsec_set.nuf(g,ct)/xsec_set.phiv(g,ct)
            xsec_set.abs(g,ct)=xsec_set.abs(g,ct)/xsec_set.phiv(g,ct)
            do gp=1, ngcd
                xsec_set.sct(gp,g,ct)=xsec_set.sct(gp,g,ct)/xsec_set.phiv(g,ct)
            enddo
            xsec_set.rmv(g,ct)=xsec_set.tot(g,ct)-xsec_set.sct(g,g,ct)
        enddo

        ! name -> xt(ct)a(iasmtype)p(pt)
2021    format ("spectrum_xt(",i<MAX(0,0)>,")A(",i<MAX(0,0)>,")P(",i<MAX(0,0)>").txt")
        write(sfile, 2021) ct, iasmtype, xsec_match(k,iasmtype).key
        open (2020, file=sfile, status="unknown")
        write(2020, '(2000(es12.5,1x))') xsec_set.phiv(:,ct)
        close(2020)
        
        
        xsec_set.chi(:,ct)=0
        
        if ( not(isfuelasm(iasmtype))) exit
        
        if (homtype .eq. HT_PTASM) then
            do g=1, ng
                xt=xsec_match(k,iasmtype).key
                xsec_set.chi(gtog(g),ct)=xsec_set.chi(gtog(g),ct)+xchi(xt,g)
            enddo
            xsec_set.chi(:,ct)=xsec_set.chi(:,ct)/sum(xsec_set.chi(:,ct))
        elseif (homtype .eq. HT_ASMHM) then
            asmpsi=asmpsi/sum(asmpsi)
            do j=1, npy
                do i=1, npx
                    ! asmtype == ct
                    xt=assembly(i,j,iasmtype)
                    do g=1, ng
                        xsec_set.chi(gtog(g),ct)=xsec_set.chi(gtog(g),ct)+xchi(xt,g)*asmpsi(i,j)
                    enddo
                enddo
            enddo
            xsec_set.chi(:,ct)=xsec_set.chi(:,ct)/sum(xsec_set.chi(:,ct))
        endif
        
    enddo 
    
endsubroutine 

function setupasm_nxt() result(ret)
    include "pspec.FI"
    integer :: ret, tasm, i, j
    logical :: findit
    findit = .false.
    iasmtype=iasmtype+1
    if (iasmtype > ntypeasm) then
        ret=0
    elseif (sum(assembly(:,:,iasmtype)) .eq. 0) then
        ! dummy
        ret=0
    else
        ret=iasmtype
    endif
    
    ! check if used in the core
    do j=1, nay
        do i=1, nax
            tasm=core(i,j)
            if (tasm.eq. iasmtype) then
                findit=.true.
                exit
            endif
        enddo
        if (findit) exit
    enddo      
    if (findit) then
        if (isfuelasm(iasmtype)) then
            call setfuelnxt()
        else
            call setreflnxt()
        endif
    else
        ret=0
    endif
endfunction

subroutine setfuelnxt()
    use define
    include "pspec.FI"
    asmgeo.ncmx=npx
    asmgeo.ncmy=npx
    allocate(asmgeo.hacx(0:asmgeo.ncmx))
    allocate(asmgeo.hacy(0:asmgeo.ncmy))
    allocate(asmgeo.cplst(asmgeo.ncmx,asmgeo.ncmy))
    asmgeo.hacx=hacx(0:asmgeo.ncmx)
    asmgeo.hacy=hacy(0:asmgeo.ncmy)
    asmgeo.nfmx=asmgeo.ncmx
    asmgeo.nfmy=asmgeo.ncmy
    asmgeo.albedo=0.0_8
    asmgeo.cplst(:,:)=assembly(:,:,iasmtype)
    asmgeo.ibndry(DW)=1
    asmgeo.ibndry(DE)=asmgeo.nfmx
    asmgeo.ibndry(DN)=asmgeo.nfmy
    asmgeo.ibndry(DS)=1
    
    
    ncmx=asmgeo.ncmx
    ncmy=asmgeo.ncmy
    nfmx=asmgeo.nfmx
    nfmy=asmgeo.nfmy
    albedo=asmgeo.albedo
    cplst=>asmgeo.cplst
    hacx=>asmgeo.hacx
    hacy=>asmgeo.hacy
    ibndry=asmgeo.ibndry
endsubroutine

subroutine setreflnxt()
    use define
    implicit none
    include "pspec.FI"
    include "xsec.FI"
!    logical :: ret
    integer :: i, j, tasm
    integer :: reflasm, mxfreq, mxcntasm
    integer :: cntnghrefl(ntypeasm)
    
    
!    ret=.true.
!    
!    if (nax < 2) then
!        ! ignore
!        ret=.false.
!        return
!    endif
    
!    reflasm = 0
!    ! search for reflector
!    ! assumint that there's only an assembly type assigned for reflector
!    do j=1, nay
!        do i=1, nax
!            tasm=core(i,j)
!            if (tasm>0) then
!                if(isfuelasm(tasm) .eq. .false.) then
!                    reflasm=tasm
!                    exit
!                endif
!            endif
!        enddo
!        if (tasm>0) exit
!    enddo        

    reflasm = iasmtype
    
    ! find most frequent asm which is neighboring with refl
    cntnghrefl=0
    do j=1, nay
        do i=1, nax
            if (isreflneigh(i,j)) then
                tasm=core(i,j)
                cntnghrefl(tasm)=cntnghrefl(tasm)+1
            endif
        enddo
    enddo
    
    mxfreq=0
    mxcntasm=0
    do tasm=1, ntypeasm
        if (cntnghrefl(tasm)>mxfreq) then
            mxfreq=cntnghrefl(tasm)
            mxcntasm=tasm
        endif
    enddo
    
    ! setting up reflective problem 
    reflgeo.ncmx=npx*2
    reflgeo.ncmy=npx
    allocate(reflgeo.hacx(0:reflgeo.ncmx))
    allocate(reflgeo.hacy(0:reflgeo.ncmy))
    allocate(reflgeo.cplst(reflgeo.ncmx,reflgeo.ncmy))
    reflgeo.hacx=coregeo.hacx(0:reflgeo.ncmx)
    reflgeo.hacy=coregeo.hacy(0:reflgeo.ncmy)
    reflgeo.nfmx=reflgeo.ncmx
    reflgeo.nfmy=reflgeo.ncmy
    reflgeo.albedo=0.0_8
    reflgeo.cplst(1:npx,:)=assembly(:,:,mxcntasm)
    reflgeo.cplst(npx+1:npx*2,:)=assembly(:,:,reflasm)
    reflgeo.ibndry(DW)=1
    reflgeo.ibndry(DE)=reflgeo.nfmx
    reflgeo.ibndry(DN)=reflgeo.nfmy
    reflgeo.ibndry(DS)=1
    
    ncmx=reflgeo.ncmx
    ncmy=reflgeo.ncmy
    nfmx=reflgeo.nfmx
    nfmy=reflgeo.nfmy
    albedo=reflgeo.albedo
    cplst=>reflgeo.cplst
    hacx=>reflgeo.hacx
    hacy=>reflgeo.hacy
    ibndry=reflgeo.ibndry
    
    
contains 
    function isreflneigh(i,j)
        integer, intent(in) :: i, j
        logical :: isreflneigh
        integer :: tasmneigh
        isreflneigh=.false.
        
        tasm=core(i,j)
        if (not(isfuelasm(tasm))) return
        
        if (i>1) then
            tasmneigh=core(i-1,j)
            if (tasmneigh .eq. reflasm) isreflneigh=.true.
        endif
        if (j>1) then
            tasmneigh=core(i,j-1)
            if (tasmneigh .eq. reflasm) isreflneigh=.true.
        endif
        if (i<nax) then
            tasmneigh=core(i+1,j)
            if (tasmneigh .eq. reflasm) isreflneigh=.true.
        endif
        if (j<nax) then
            tasmneigh=core(i,j+1)
            if (tasmneigh .eq. reflasm) isreflneigh=.true.
        endif
        
    endfunction    
endsubroutine


subroutine restore_core
    implicit none
    include "pspec.FI"
    ncmx=coregeo.ncmx
    ncmy=coregeo.ncmy
    nfmx=coregeo.nfmx
    nfmy=coregeo.nfmy
    albedo=coregeo.albedo
    cplst=>coregeo.cplst
    hacx=>coregeo.hacx
    hacy=>coregeo.hacy
    ibndry=coregeo.ibndry
endsubroutine

subroutine findFASM
    implicit none
    include "pspec.FI"
    include "xsec.FI"  
    integer :: tpin, tcomp, tasm
    integer :: i, j

    ! fuel pin search
    do tpin=1, ntypepin
        tcomp = typepin(tpin)%rc(typepin(tpin)%nr)
        if (sum(Xnuf(tcomp,:)) .ne. 0) isfuelpin(tpin)=.true.
    enddo
    
    ! fuel assembly search
    do tasm=1, ntypeasm
        do j=1, npy
            do i=1, npx
                if (isfuelpin(assembly(i,j,tasm))) then
                    isfuelasm(tasm) = .true.
                    exit
                endif
            enddo
        enddo
    enddo
endsubroutine


subroutine pprec(psi, ia, ja)
    include "pspec.FI"
    integer :: ia, ja, i, j
    real(8) :: psisum, psisumafter, fmod
    real(8), dimension(npx, npy) :: psi
    real(8), pointer, dimension(:,:) :: asmpsi
    integer :: asmtype
    asmtype=core(ia,ja)
    if (isfuelasm(asmtype)) then
        asmpsi=>shapef(:,:,asmtype)
        psisum=sum(psi(1:npx, 1:npy))
        do j=1, npy
            do i=1, npx
                psi(i,j)=psi(i,j)*asmpsi(i,j)
            enddo
        enddo
        psisumafter=sum(psi(1:npx, 1:npy))
        fmod=psisum/psisumafter
        
        do j=1, npy
            do i=1, npx
                psi(i,j)=psi(i,j)*fmod
            enddo
        enddo
    endif
endsubroutine
ENDMODULE asmcal